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Abstract. We derive the analytical solutions of the bound timelike geodesic orbits in 
Kerr spacetime. The analytical solutions are expressed in terms of the elliptic integrals 
using Mino time A as the independent variable. Mino time decouples the radial and 
polar motion of a particle and hence leads to forms more useful to estimate three 
fundamental frequencies, radial, polar and azimuthal motion, for the bound timelike 
geodesies in Kerr spacetime. This paper gives the first derivation of the analytical 
expressions of the fundamental frequencies. This paper also gives the first derivation 
of the analytical expressions of all coordinates for the bound timelike geodesies using 
Mino time. These analytical expressions should be useful not only to investigate 
physical properties of Kerr geodesies but more importantly to applications related 
to the estimation of gravitational waves from the extreme mass ratio inspirals. 
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1. Introduction 

The Kerr black hole has been well studied since the discovery of the Kerr solution. It 
is an important topic not only in mathematical problems of general theory of relativity, 
but also for applications in astrophysics. Currently, there are many candidates for black 
holes in the universe and they have a wide range of mass scales ranging from stellar 
mass scales to galactic nuclei mass scales pp. 

One of the ways to investigate the properties of a Kerr black hole spacetime is to 
study geodesic motion in this background. Detailed works on the geodesic motion in 
black hole spacetimes are summarized in Chandrasekhar [2]. In the weak field regime, 
at large distances from the black hole, the orbits of a particle are almost the same as 
that in Newtonian gravity. In the strong field regime, however, the orbits become more 
complicated and it is difficult to compare the orbits with that in Newtonian gravity. 
For the case of bound geodesies, this can be explained by mismatches between the 
fundamental frequencies of radial, Q r , polar, Qg and azimuthal-motion, Q^. For example, 
Qc/, — Qg shows the precession of the orbital plane and — Q r shows the precession of 
the orbital ellipse. Differences between the fundamental frequencies become larger as 
the particle goes into the strong gravity region around black hole horizon or separatrix, 
which is the boundary between stable and unstable orbits. These relativistic effects 
have been studied for some cases and some examples of extreme phenomena are found 
as follows. 

Wilkins [3] derived the analytical expressions for the ratio of the azimuthal 
frequency and the polar frequency, fl^/Qg, when a particle moves on both circular 
and non-equatorial orbits around the extreme Kerr black hole. He then showed that the 
ratio becomes larger as the particle approaches the horizon and found that the particle 
traces out a helix-like orbit on a sphere around the black hole. He also pointed out 
that there exist horizon- skimming orbits which have the same radius as the horizon. 
Horizon-skimming orbits are also studied by numerical calculations including the effects 
of the emission of gravitational waves from a particle for circular and non-equatorial 
orbits [1] and for generic orbits [5] around near-extremal Kerr black holes. Glampedakis 
and Kennefick [6] numerically investigated the ratio of the azimuthal frequency and the 
radial frequency, f2^/f2 r , when a particle moves both on eccentric and equatorial orbits 
around the Kerr black hole. They found that the ratio becomes larger as the particle 
approaches the separatrix and the particle traces out a quasi-circular orbit around the 
periapsis before going back to the apoapsis. These orbits are called zoom-whirl orbits. 

The above results show that the fundamental frequencies play an important role in 
understanding bound geodesic orbits. However, the coupling of the r and 6 motions in 
the geodesic equation has prevented one from deriving the fundamental frequencies, Q r , 
Qg and f^, for general bound geodesic orbits until recently. Using the elegant Hamilton- 
Jacobi formalism, Schmidt [7] derived the fundamental frequencies without discussing 
the coupling of the r and 9 motions. Although his results show that we can expand 
an arbitrary function of the particle's orbit in a Fourier series, we can not estimate the 
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Fourier components because of the coupling of the r and ^-motion. Mino [8] showed that 
we can separate r and ^-motion if we use new time parameter A and derived the integral 
forms of the periods of both r and ^-motion with respect to A, which is called Mino 
time. Combining Schmidt's method with Mino time, Drasco and Hughes [9] derived the 
fundamental frequencies and showed how the Fourier components of arbitrary functions 
of orbits with respect to Mino time can be computed because of the decoupling of both 
r and 6 motions. They also showed how from these results using Mino time, the Fourier 
components with respect to coordinate time can also be derived. Thanks to these 
results, one can compute gravitational waves from binary systems in which a stellar 
mass compact star is moving on a general bound geodesic orbit around a supermassive 
black hole, the so-called extreme mass ratio inspirals(EMRIs) [10]. Gravitational waves 
from EMRIs are one of the main targets for space-based Laser Interferometer Space 
Antenna (LISA) p]. 

In this paper, we derive analytical expressions for bound timelike geodesic orbits 
in Kerr spacetime using Mino time as the independent variable. Despite a lot of works 
on geodesic motion [2], the analytical expressions of null or timelike geodesies in Kerr 
spacetime are still important subjects. Fast and accurate computation of null geodesies 
in Kerr spacetime is required to study radiation which pass near black holes in accretion 
systems such as active galactic nuclei and X-ray binaries (see, for example, [121 13] and 
references therein). Fast and accurate computation of timelike geodesies is also required 
to study gravitational waves from EMRIs and construct efficient templates for LISA 
data analysis. Rauch and Blandford gave tables which reduce the some integral forms 
of null geodesies to Legendre elliptic integrals p3] using the radial coordinate as the 
independent variable [12]. They did not give the complete tables which reduce all the 
integral forms to the elliptic integrals because it was easier and faster to compute both 
t and (f) coordinates numerically when they studied the optical structure of the primary 
caustic around Kerr black hole. Using Carlson elliptic integrals pj)] to calculate all 
coordinates of null geodesies, however, Dexter and Agol showed that they can compute 
null geodesies more efficiently than numerical integration method |13j . Although they 
did not show analytical expressions of all coordinates of null geodesies since there are so 
many cases to be considered, they opened their numerical code to compute null geodesies 
semi-analytically in Kerr spacetime. In this paper, we show that we can easily derive the 
analytical expressions of bound timelike geodesies in terms of Legendre elliptic integrals 
if we properly transform the r and 9 variables. This is the first time that the analytical 
expressions of fundamental frequencies are derived. This is also the first time that the 
analytical expressions of all geodesic coordinates are derived using Mino time as the 
independent variable. These analytical expressions of bound timelike geodesic orbits 
with respect to Mino time are simpler than that in [13] for null geodesies and should be 
useful to investigate gravitational waves from EMRIs. The analytical solutions should 
also be helpful for investigations of bound geodesies in Kerr spacetime. 

This paper is organized as follows. In section [2j we review Kerr geodesies using 
observer time. We then discuss Kerr geodesies in Mino time and derive the analytical 
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expressions of the fundamental frequencies of bound geodesies in section [3J In section HJ 
we derive the analytical expressions for bound geodesic orbits. We conclude with a 
brief summary in section [5j In the Appendices, we discuss technical details of the 
implementation required to obtain the results in this paper. Throughout this paper, we 
use units with G — c — 1. 




2. Geodesic Orbits in Kerr Spacetime 

The geodesic equations that describe a particle's orbits in Kerr spacetime are given by 

2 

= R(r), 

£ — — — ) = 6(cosi 
dt 

S3- = T r (r) + T e (cos6) + aC z , 
dr 

E-^ = $ r (r) + $ e (cos#) - aS. (1) 
dr 

The functions R(r), 6(cos#), T r (r), T^cos^), $ r (r) and ^(cos^) are defined by 

R(r) = [P(r)] 2 - A[r 2 + (a£ - C z f + C], 

0(cosfl) = C-(C + a 2 (l - £ 2 ) + £ 2 ) cos 2 + a 2 (l - £ 2 ) cos 4 0, 

2 2 

T r (r) = r ^ P(r), T e (cos#) = -a 2 £(l - cos 2 #), 

$ r (r) = ^P(r), $ e (cos^) = - £ % 
A 1 — cos z a 

with P(r) = £(r 2 + a 2 ) - a£ z , S = r 2 + a 2 cos 2 # and A = r 2 - 2Mr + a 2 . Here M and a 

are the mass and the angular momentum of the black hole, respectively. There are three 

constants of motion, S, C z and C, which are the energy, the z-component of the angular 

momentum and the Carter constant per unit mass, respectively. Using reasonable initial 

conditions for the particle's orbit, we can derive the orbits using the proper time of the 

particle, r, by numerical integration. Dividing dr/dr, dcos^/dr and d0/dr by dt/dr, 

we can also derive the orbits with coordinate time t by numerical integration. When the 

orbits are bound to black hole, however, we have to take care of the turning points in 

the radial and the polar motion where the signs of dr/dr and dcos#/dr change. These 

turning points correspond to periapsis and apoapsis for the radial motion, and 6* min and 

7T — #min for the polar motion, where # m i n is the minimum value of 9. We need smaller 

stepsizes to resolve the derivatives around turning points. We can avoid this problem 

by introducing new variables for the radial and polar motion, r = pM/ (1 + e cos ip) and 

cos# = cos 6* m i n cos x, where p is semilatus rectum and e is eccentricity [9]. Using these 

new variables, ip and x, we can estimate the orbits accurately. 

There exists three fundamental frequencies, Q r , Q$ and Q^, for bound Kerr 

geodesies. However, it is difficult to estimate the fundamental frequencies using ([1]) 

because of the coupling of the r and ^-motions. For instance, we immediately face 
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a difficulty when we estimate fl r using dr/dt = (dr/dr)(dt/dr) 1 because r and 9 
asynchronously pass their turning points. 



3. Fundamental Frequencies of bound geodesies 

We now proceed to derive the analytical expressions for the fundamental frequencies 
of bound geodesic orbits using Mino time. In section 13.11 we briefly describe the Kerr 
geodesies in Mino time and then show how to derive the analytical expressions for the 
fundamental frequencies in section I3~2l and section [3731 In section l374"t we will check the 
analytical expressions by comparing them with earlier literature. 



3.1. Geodesies in Mino time 



Using Mino time, A = / dr/S, the geodesic equations become 

2 

I =R(r), 

2 

| =e(cos#), 



( dr 
ydx) " 

'd cos s 



zi 



\ dX ) 
dt 

— = T v (r)+T e (cos6) + a£ 
dA 

^ = $ r (r) + $ e (cos#)-a£. (2) 

It should be noted that, in fl2j), dr/dX depends only on r and dcos^/dA depends 

only on cos 6. Thus the equations for radial and polar motion are decoupled. For the 

bound orbits, r(A) and cos#(A) become periodic functions which are independent of 

each other. The fundamental periods for the radial and polar motion, A r and Ag, with 

respect to A are given by 

f™ dr . rc°s0mm dcos6> 

A r = 2 / -==, A, = 4 / (3) 
Jr min y/R( r ) Jo ^e(cosfl) 

where 

PM pM 7T 

inc + (sgn£J 9 min = -. (4) 



■ Him -| i ■ ijjlcia. -f j nit, 1 \ O */ iiiiii r% 

1 + e 1 — e 2 

Here r min and r max are the periapsis and apoapsis for the radial motion respectively, and 
8 inc the inclination angle from the equatorial plane of black hole. Of course, (S,C Z ,C) 
are described by these orbital parameters (p, e, 9i nc ) and given in [7J [9]. The angular 
frequencies of the radial and the polar motion then become 

We also note that both dt/dX and d(f>/dX in (j2J) are the sum of a function of r and 
a function of cos 9. Then each equations are integrated as 

t{X) = TX + &\X) + i w (A), 0(A) = T^A + (r) (A) + W (A), (6) 
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where T and are the frequencies of coordinate time t and with respect to A 
respectively, which are given by 

r = T t ( r ) + T t (6) + aC z , = T^(r) + T^e) — a£, 

T>) = (T r (r)) A , T tW = (T e {cos6)) x , 

T#r) = (9 r (r))x> TC> = (7) 

where (• ■ -) A = liniAA^oo(2AA)~ 1 j A ^ x dX - ■ ■ represents infinite time average with respect 
to A, and and 0^/^ satisfy 

dt (r) dt^ 

—— = T r (r) - T tW , -— = T e (cos0) - T tW , 

^- = $ r (r) - T>), = $ e (cos#) - T>. (8) 

Equation ([6]) shows that both t(X) and 0(A) consist of two distinct parts. The 
first term represents an accumulation over A-time and the last two terms represent 
oscillations around it with periods 2n/T r and 2-k/Tq. We note that the frequencies 
with respect to A are related to the frequencies with distant observer time as [9] 

r-^*, i-^, 

&r = -yf , = -jr, &><f> = -yr- (9) 

In the following subsections, section 13.21 and section I3.3[ we discuss the analytical 
expressions for these frequencies. And we discuss the analytical expressions of the orbits, 
r(A), cos6'(A), t(A) and 0(A), in section HI 

3.2. Frequencies of r andO-motion 

In this subsection, we derive the analytical expressions for the frequencies of r and 9- 
motion, T r and Tg, using ([3]). As explained in section^ R(r) and 6(cos#) become zero 
when r and cos 6* go through the turning points, r min , r max and icos^min, respectively. 
Thus we usually transform the variables, r and cos 8, to avoid divergences in the 
numerical calculation. However we know that can be expressed in terms of the 
elliptic integrals since both R(r) and 0(cos#) are fourth order polynomials [16]. It is 
useful if we know the four zero points of both R(r) and 0(cos#) to express ([3]) in terms 
of the elliptic integrals. We rewrite R(r) and 0(cos#) as [9] 

R(r) = (1 — £ 2 ){ri — r)(r — r 2 )(r — r 3 )(r — r 4 ), 

6(cos#) = C 2 z e (z_ - cos 2 6)(z + -cos 2 6), (10) 

where 



pM pM (A + B) + J(A + B) 2 -AAB AB 

ri = —' r2 = it? r3 = 2 — ; ' r4 = — ' 

A + B-JM;-^), ^= (1 _g rira > (ID 

and where e = a 2 {l — £ 2 )/£ 2 z , z_ = cos 2 9 min and z + = C/(C 2 z e z_). We note that 
two zero points, r\ and r 2 , of Rir) are apoapsis and periapis respectively and two zero 
points, 2_ and —z_, of 9(cos6 l ) are # min and n — 9 min respectively. These zero points 
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correspond to turning points, defined in of radial and polar motion. But the other 
two zero points of both R(r) and 0(cos6 l ), r 3 , r 4 and iz + , do not correspond to turning 
points of radial and polar motion. 

Using ({TO]) , we can express ([3]) in terms of the elliptic integrals as 
dr' 2 

=F(arcsiny r , k r ), 



■<*2 yjR(r') yj(l -S 2 )^ -r 3 )(r 2 -r 4 ) 

rcos9 dcOS0' 1 . , , , . 

/ / = 7 — ^^^(arcsinye, fc e ), (12) 

J o J Q (cos 8') £zy/eo^T 



where 



r 1 -r 3 r-r 2 r 1 - r 2 r 3 - r 4 

V Ti — r 2 r — r 3 V r i — r 3 r 2 — r 4 

Ve = ^=, k e = \ — , (13) 



and .F(<^, fc) is the incomplete elliptic integral of the first kind defined by 

F(<p,k)= -. = V (14) 

Jo Vl-A; 2 sin 2 # io /(i-^)^.^) 



In the following, we describe both the elliptic integrals and the elliptic functions using 
the notation in [17J. The orbital frequencies of radial and polar motion with respect to 
A are then given by 



2K(k r ) ' 9 2K{k e ) ' 1 ' 

Here K(k) is the complete elliptic integral of the first kind defined by K(k) = F(tt/2, k). 
We note that the analysis of ^-motion here is similar to that of Drasco and Hughes [9]. 
Though they used a different transformation of cos#, our expression for Tg in this 
subsection agrees with their final result. 

3.3. Frequencies oft and <p -motion 

In this subsection, we derive the analytical expressions for the frequencies of t and 4>- 
motion, T and T^, using (17j) . Since the r and ^-motion decouple in Mino time, we can 
rewrite the infinite time average in ([7j) as an average over an orbital period, A r or Ag, 

as 

2 rn T r (r) 4 f V*= Tg(cos9) J 

T t (»-) = — / dr, T t ( 9 ) = — / =dcosfl, 



A r^ 2 J_R( r ) ' ' A Jo Je(cos^) 



2 rn $ r r 4 $ e (cos0) 

T <^w = X" / r— dr, Tt>) = — / dcosfl. 16 

A r^2 ./R( r ) Ae^o (cos 0) 



It is straightforward to express Y t (o) and Tj,(«) in terms of the elliptic integrals if we use 
ye in (El). 

2a 2 ST e 



T 



t w 



7T 



'l-Z + )K(kg)+Z + E(^,kg) 
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T^(fl) = — — — n(-, fo), (17) 

where -E(y2, fc) is the incomplete elliptic integral of the second kind and H(ip, c, k) is the 

incomplete elliptic integral of the third kind defined by 

r<e / 

,2, 



E(tp,k) = / Jl- k 2 sin 2 9d9, 
Jo 

rf d9 

n(v,c,* =/ 7: . 2m 71 /2 . 2 =- (is) 

Jo (1 + csm 2 9)Vl - k 2 sm 2 9 
Note that, E(n/2, k) is the complete elliptic integral of the second kind and II(7r/2, c, fc) 
is the complete elliptic integral of the third kind. In the followings, we describe E(ir/2, k) 
as E(k), and H(tt/2, c, k) as U(c, k). 

On the other hand, we have to rewrite T r {r) and $ r (r) in order to express T t ( r ) 
and T^w in terms of the elliptic integrals. Performing partial fraction decomposition, 
we decompose T r {r) and $ r (r) as follows. 

T r (r) = Er* + 2M£r + -™- {(^ - ^.W - 2Jfq»£ ' 

+ (a 2 + 4M 2 )^-a£ 2 (|a| ^ M), 

^2 2M(AM 2 £-aC z ) 2M 2 (2M 2 £ - aC z ) 

= Sr 2 + 2M£r + - 



r — M (r - M) 2 

+ (a 2 + 4M 2 )f-a£ z (|a|=M), 



a 



{ 2M£r + - aC 



$ r (r) = i - - (+ <- -H + a£ (\a\ ? M), 

r + — r_ [ r — r + J 

2Ma£ a(2M 2 £-aC z ) „ ,. . 



with r± = M ± V 'M 2 — a 2 . From (fl3l) . we find r = t 3 + (r 2 — r 3 )/(l — h r y 2 ), where 
= ( r i — r 2) I iji — r z)- Then it is straightforward to compute J r r 2 r'dr' / \jR(r'), which 
is composed in T t (V), as 

r r' ^ _ 2 [r 3 F(arcsiny r , fc r ) + (r 2 - r 3 )n(arcsiny r , -fe r , k r )} _ 

Jr* Jnjy) r ^(l-^Xn-r,)^-^) 



We show the results of the other terms of T r (r')dr' / J R(r') and JT § r (r')dr' / J R(r') 



in 



Appendix A Using results quoted there, we derive T and as 



r = am 2 £ + 2a2£z ^_ [K{ke) _ s^)] + 2T - 

tt£ 2V ^7 tt^CI -£ 2 )(ri -r 3 )(r 2 - r 4 ) 

x | g t( r 3( r i + r 2 + r 3) - r l r 2 )K(k r ) 
+ (r 2 - r z )(r 1 + r 2 + r 3 + r 4 )n(-/i r , fc r )] 
+(ri - r 3 )(r 2 - r 4 )£(£; r ) + 2M£ [r 3 K(k r ) + (r 2 - r 3 )n(-/i r , fc, 
2M 



+ 



r — r_ 



r 3 - r + \ r 2 - r + 

-(+--)]}, 

t, = -^H-z-M + , 2aTr 

^V^T vr(r + - r_) v /(l - £ 2 )( ri - r 3 )(r 2 - r 4 ) 
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r 3 



K(k r ) - 



r2 ~ r 3 



II(— h + , k r 



- + 



(21) 



r 2 - r+ 

where h± = (ri — r 2 )(r 3 — r±)/[(ri — r 3 )(r 2 — r±)]. Combining Qj, f|T5l) and fl2TT) . we can 
derive the orbital frequencies with respect to observer time, Q r , Qg and fl^. It should 
be noted that ( 12TI) is not valid for the case \a\ — M since there exists divergent terms in 
l/(r + — r_). We show the expressions of T and for the case \a\ = M in|Appendix B 



3.4- Consistency check of the fundamental frequencies 

We can compare the expressions for the fundamental frequencies, fl r , fig and fl^, in this 
section with that in earlier literature for some limiting cases. We can compare fl^ for 
the case e = and 9 inc = with [18], fl^/flj. in the case a = 0, e ^ and 9 inc = with 
[T9] , and both fl^ and fig in the case e = and 6 inc ^ with |20j . It is a good check on 
our results that our more general analytical expressions for the fundamental frequencies 
are consistent with earlier work for the limiting cases. 

We can check our results for more general cases. In table [Hand table EJ we compare 
our results with [21] in which the analytical expressions of the orbital frequencies are 
derived in terms of both post-Newtonian and small eccentricity expansions up through 
0(v 5 , e 2 ), where v 2 = M/p. In table [H and table EJ we check our results for the cases 
a ytz M and a = M respectively. We find that relative errors are always less than 10 -4 
when we compare the results for the cases p = 100M and 0.01 < e < 0.09. Since the 
fundamental frequencies in [21] are derived up through 0(v 5 , e 2 ), these relative errors, 
less than 10~ 4 , show the consistency of our results with [21]. In table El we compare 
our results with numerical integration method for the eccentric and inclined orbits such 
that p = 6M, e = 0.7, #i nc = 20° and a = 0.9M or a = M. In numerical integration 
method, we use the trapezium rule to compute the fundamental frequencies. We can 
compute very accurately if we use the trapezium rule for the numerical integration of 
a periodic function. Then we find that the analytical expressions of the fundamental 
frequencies in this section agree with the results of numerical integration method. The 
relative errors are less than 10~ 15 in double precision calculation. These facts show that 
the analytical expressions in this section are correct in the cases of generic bound orbits. 



4. Analytical solutions of bound geodesies 

In this section, we derive the analytical expressions for bound geodesic orbits, r(A), 
cos#(A), t(X) and 0(A), in terms of the elliptic integrals. Since we have already derived 
the orbital frequencies in terms of the complete elliptic integrals in section [31 we can 
derive the orbits if we replace the complete elliptic integrals with the incomplete elliptic 
integrals. However, we have to take account of the initial values of both r and 9 and 
the signs of both dr/dA and dcos#/dA at given A-time. In the following subsections, 
we derive the radial solutions, r(A), t^(A) and 0^(A), in section fl~Tl and the polar 
solutions, cos#(A), t^ e \\) and 0^(A), in section I4T21 Finally, we check the consistency 
of our analytical results by comparison with the results of earlier literature in section l4~3l 



Analytical solutions of bound timelike geodesic orbits in Kerr spacetime 



10 



Table 1. Comparison of the orbital frequencies, Q r , fig and fl^, derived using 
analytical expressions in this work and the analytical post-Newtonian expressions for 
orbits which are slightly eccentric but greatly inclined [3T] in the case of a = 0.9M and 
p = 100M. Our results are consistent with post-Newtonian results. Relative errors of 
the orbital frequencies are always less than 10~ 4 . 
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Table 2. Comparison of the orbital frequencies, fl r , fie and fi^, derived using 
analytical expressions in this work and the analytical post-Newtonian expressions for 
orbits which are slightly eccentric but greatly inclined [3T] in the case of a — M and 
p = 100M. Our results are consistent with post-Newtonian results. Relative errors of 
the orbital frequencies are always less than 10 -4 . 



e 


Sine 






ork 








ork 




nj his 

V 




ork 




Q p ost - 


Newto 




Post — Newto 




n Post- 


Ne 






0.01 
0.01 
0.01 


20° 
45° 
70° 


9.72213 
9.71546 
9.70462 


X 
X 
X 


10" 

io- 

10" 


4 
4 
4 


9.97159 
9.97769 
9.98777 


X 
X 
X 


10" 
10" 
10" 


4 
4 
4 


9.99017 
9.99663 
1.00073 


X 
X 
X 


10" 
10" 
10" 


4 
4 
3 


9.72250 
9.71575 
9.70484 


X 
X 
X 


10" 

io- 
io- 


4 
4 
4 


9.97157 X 
9.97770 X 
9.98780 X 


io- 

10" 

io- 


4 
4 
4 


9.99017 
9.99664 
1.00073 


X 
X 
X 


10" 
10" 
10" 


4 
4 
3 


0.05 
0.05 
0.05 


20° 
45° 
70° 


9.68780 
9.68117 
9.67039 


X 
X 
X 


10" 
10" 

io- 


4 
4 
4 


9.93637 
9.94247 
9.95254 


X 
X 
X 


10" 
10" 
10" 


4 
4 
4 


9.95490 
9.96134 
9.97196 


X 
X 
X 


10" 
10" 
10" 


4 
4 
4 


9.68827 
9.68157 
9.67072 


X 
X 
X 


10" 

io- 
io- 


4 
4 
4 


9.93647 X 
9.94259 X 
9.95268 X 


io- 

10" 

io- 


4 
4 
4 


9.95501 
9.96147 
9.97210 


X 
X 
X 


10" 
10" 
10" 


4 
4 
4 


0.09 
0.09 
0.09 


20° 
45° 
70° 


9.60784 
9.60130 
9.59067 


X 
X 
X 


io- 
io- 

10 _ 


4 
4 
4 


9.85436 
9.86045 
9.87050 


X 
X 
X 


10" 
10" 
10" 


4 
4 
4 


9.87273 
9.87916 
9.88975 


X 
X 
X 


10" 
10" 
10" 


4 
4 
4 


9.60935 
9.60274 
9.59204 


X 
X 
X 


io- 

10" 
10" 


4 
4 
4 


9.85553 X 
9.86164 x 
9.87171 x 


io- 
io- 

10" 


4 
4 
4 


9.87391 
9.88036 
9.89097 


X 
X 
X 


10" 
10" 
10" 


4 
4 
4 



Table 3. Comparison of the orbital frequencies, Q r , fig and fl^, derived using 
analytical expressions in this work and numerical integration method in the case of 
p = 6M, e = 0.7, 9 U1C = 20° and a — 0.9M or a — M. Our results are consistent with 
numerical integration method. Relative errors of the orbital frequencies agree with the 
accuracy of double precision calculation. 



a/M 




This work 




Numerical integration 




Absolute value of relative error 


0.9 


n r 


1.8928532285101992 x 10" 


■2 


1.8928532285101982 x 10" 


■2 


5.6 X 10" 


16 


0.9 


n e 


2.7299110395017517 x 10" 


2 


2.7299110395017506 x 10" 


-2 


4.1 X 10" 


16 


0.9 


n,f, 


3.0550463796964692 x 10" 


2 


3.0550463796964682 x 10" 


-2 


3.6 x 10- 


16 


1 


n r 


1.9343466898960462 x 10" 


■2 


1.9343466898960444 x 10" 


■2 


7.8 x 10" 


16 


1 


Oh 


2.6337035996626332 x 10" 


2 


2.6337035996626321 x 10" 


-2 


4.3 x 10" 


16 


1 


0<f, 


2.9662029663040452 x 10" 


2 


2.9662029663040452 x 10" 


-2 


4.0 x 10" 


17 
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4-1. Radial solution :r(A) , t^(X) and (j)^ r \\) 

Solving ([12]) and (ED, we obtain A(r), &\\) and (r) (A) as 

A(r) = 



dr' 



R(r') 



{r) (A) 



•(A) T ( r ') - T 



t(0 



i2« 

r(A) $ r ( r ') - L(r) 



dr', 



-R(r') 



dr'. 



(22) 



We derive r(A) inverting A(r). Since the period of r-motion with respect to A is 
A r = 27r/T r , we map A to AM as A^ = A - 27r[T r A/2vr]/T r , where [• • •] is the floor 
function, in the following subsections. In order to investigate the integrations in (1221) 
properly, we have to take account of r(A = 0) and the sign of dr(A)/dA. There exist two 
cases depending on whether the initial value is dr(0)/dA > or dr(0)/dA < 0. In the 
following subsections, we consider the two cases separately. We note that the expressions 
of both t^ r \\) and (f)^ T \X) in the following subsections are valid when \a\ 7^ M. We show 
t^ r \X) and (j)( r \X) when \a\ — M in Appendix B 



4-1.1. dr(0)/dA > case In this subsection, we consider the case that the initial value 
of r(A) satisfy dr(0)/dA > 0. We set r(A = 0) = Tq in this subsection. Then A(r) in 
( |22l) can be expressed as 

rr cl r ' 



R{r') 

(i) 



r 2 Jr 2 



dr' 



r : r. 



(i) 



R{r') 



n, 



'1 



^ / +2 

r 2 J t°2 J r 2 



J 1 ) 



dr' 



R{r') 



+2 

T2 JT2 J T2 



AD 




dr' 



R{r') 



Thus we find the solution as 

f Aj r, (r)-AW 
A«(r) = | _A^(r)+A r -AW 
{ Af ) (r)+A,-AW 

where 



r : r 



r : ri — > r 2 , 



(i) 

r : r2 — >• Tq . 



^ 'i) 



and A( 1 ) = Aj ) r) (r( 1) ). 



1 



V 7 ! - £ 2 J (n - r 3 ) (r 2 - r 4 ) 



r : r x — > r 2 , 
(i) 

r : r 2 — > rg , 



F (arcsin?/ r , fc r ) , 



(23) 



(24) 



(25) 
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Inverting (1241) . we derive r(A) as 

m _ r 3 (n - r 2 )sn 2 (M r (A), fc r ) - r 2 (ri - r 3 ) 
n) ~ (r 1 -r 2 )sn2(n r (A),A; r )-(r 1 -r3) ' 1 } 

where sn(ii, k) is Jacobi's elliptic function which is defined as the inverse function of the 
incomplete elliptic integrals, u = F(<p, k), and 

[ 2K(k r )(\W + A«)/A r (0< AW <A r /2-A«), 

u r (X) = I 2K(kr)(-\W + A r - A«)/A r (A r /2 - A« < A< r > < A r - A«), (27) 

[ 2K(k r )(X^ - A r + A«)/A r (A r - A« < A« < A r ). 

Combining the results of f[2"5|) with the results of Appendix A, we can derive t^) 
and (f) {r) in ([22]) as 

2 



(r) 



^(1 - ^ 2 )(r x - r 3 )(r 2 - r 4 ) 
£ r 

x 



(r- 2 - r 3 )(ri + r 2 + r 3 + r 4 )II r (V'r, -Ar, K) 
+ (n - r 3 )(r 2 - r 4 )E r (ip r , h r , k r ) 

+2M£(r 2 - r 3 )n r (Vv, -fc r , fc r ) 

(4M 2 £ - aC z )r + - 2Ma 2 £ r 2 - r 3 ~ 



(r) 



2M 

r + — r_ 

-(+-")]}, 
2a 



r 3 - r+ 



r 2 - r + 



(2Mgr + -a£,)(r 2 -r 3 ) ~ 
(r 3 - r + )(r 2 - r+) 



(r+ - r_)y / (l -£ 2 )(r! -r 3 )(r 2 - r 4 ) 

-(+--)], (2f 
where ^ r = arcsin[sn(w r , fc r )], E r (ip r , c, fc r ) = E r (ijj r , c, k r ) — TrX r E(k r ), tl r (ip r} c,k r ) = 
T\. r {ip r , c, fc r ) — YrA r) II(c, fc r ) and 



sinz/vv (1 — sin ib r )(l — fc 2 sin ib r ) 

£ 0) (i, c, fc r = E(ip r , k T ) + — . 2 , - i r — , 

sin ip r — c 1 

E r {lP r , C, k r ) = £?^(Vr, C, fc r ) " ^ 0) (<M0), C, fc r ) 

for < A (r) < A r /2-A«, 
= -£( 0) (</v,c,fc r ) + 2£(fc r ) -E( 0) (^ r (0),c,fc r ) 

for A r /2 - A« < A (r) < A r - A«, 
= £< 0) (Vr, c, Av) + 2£(fc r ) - ^ 0) (</v(0), c, fc r ) 

for A r - A« < A (r) < A r , 

rir(Vv> c, fc r ) = n(y> r , c, fc r ) — n(^ r (o), c, fc r ) 

for < A (r) < A r /2 - A r 1} , 

= - n(Vv> c, fc r ) + 2n(c, fc r ) - n(-0r(o), C, fc r ) 

for Ar/2 - A« < A {r) < A r - A r 1} , 

= II(V>r, C, fc r ) + 2II(c, fc r ) - n(^r(0), C, fcr) 

for A r - A« < A (r) < A r . (29) 
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4-1.2. dr(0)/dA < case In this subsection, we consider the case that the initial value 

(2) 

of r(A) satisfy dr(0)/dA < 0. We set r(A = 0) = rjj in this subsection. Then A(r) in 
( |22l) can be expressed as 

rr dr' 



A (r) (r) 







R(r') 



r 2 JT2 

,(2) 


+ " 

r 2 J r 2 



dr' 



R(r') 



dr' 



R(r') 



- 1 +2 7 + 



,(2) 




dr' 



-R(r') 



(2) 

r : rg — > r 2 , 



r : r 2 — > r x , 



(2) 

r : ?"i — > 7q . 



(30) 



Thus we find the solution as 

-Af } (r)+A( 2 ) 
\ {r \r) = { A< r) (r)+Af 

-A^(r)+A r + A( 2 ) 



(2) 

r : To ' — > r 2 , 
r : r 2 — > n, 

(2) 

r : T\ — ► Tq , 



(31) 



where = X^\r^). 

Then we obtain r(A) in the same form in ( l26i) inverting ( 13TI) . but it r (A) in ( 1261) is 
modified as 

[ 2K(A; r )(-AM+A( 2 ))/A r (0 < A« < A( 2 )), 

u r (X) = I 2K(kr)(\V - A( 2 ))/A r (A( 2 ) < A« < A r /2 + A( 2 )), (32) 

[ 2K(k r )(-\^ + A r + A(, 2 ))/A r (A r /2 + A,[ 2 ) < A (r) < A r ). 

Combining the results of (|5U|) with the results of Appendix A, we can derive t^ 
and (fyrt in the same form in (T28]) . but E r (ip r , c, fc r ) and n r ( , r , c, fc r ) are modified as 

E r (ll> r , C, fc r ) = - E^fa, C, fcr) + ^(Vr(0), C, fc r ) 

for 0<A (r) <A< 2) , 

= £{°>(W,C, fc r )+^r (0) (^(0),C, fc P ) 

for A( 2 ) < A (r) < A r /2 + A^, 

= - ^° } (VV, C, Av) + 2E(K) + £< 0) (Vv(0), C, kr) 

for A r /2 + A< 2) < A (r) < A r , 

n r (^ r , c, fc r ) = — n^j., c, fc r ) + n('0 r .(o), c, fc r ) 

for < A (r) < A< 2) , 

= n(Vv,c, fc r ) + n(v» r (o),c, fc r ) 

for A^ 2) < A (r) < A r /2 + A^, 2) , 
= - U(ip r , c, k r ) + 2IT(c, kr) + n(Vv(0), c, fc r ) 

for A r /2 + A( 2) < A (r) < A r . (33) 
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4.2. Polar Solution: cosO (A), *W(A) and 0^(A) 
Solving (DGD and ©, we obtain A(cosfl), t {6) {\) and as 
X(9) - 



cos e d cos 9' 



W (A) 



e(cos^) 

cosd T e (cos9') - T t( 9) 



e(cos^') 

cos0 $J C os9') - T„ 



dcos^', 



-dcos^'. 



(34) 



^(cosfl') 

We derive cos6 l (A) inverting X(9). Since the period of ^-motion with respect to A is 
A d = 2n/T e , we map A to A^ as X^ = X — 2n[T e X/2ii}/Y e in the following subsections. 
In order to investigate the integrations in (134"|) properly, we have to take account of 
cos#(A = 0) and the sign of dcos#(A)/dA. There exist two cases depending on whether 
the initial value is dcos#(0)/dA > or dcos#(0)/dA < 0. In the following subsections, 
we consider the two cases separately. 



4-2.1. dcos9(0)/dX > case In this subsection, we consider the case that the initial 
value of cos 0(A) satisfies dcos#(0)/dA > 0. We set 9{X = 0) = 9 { 1] in this subsection. 
Then X(9) in ( 1341) can be expressed as 
rcose dcos#' 

JonoP ,/e(cos^)' 



cos 8 





cos 8 



+2 



-4 



g(l) 



COS min 



COS min 



d cos 9' 

e(cos^) 



9 : 9q^ > 9 m i n , 



cos 8. 



(1) 



d cos 9' 



10 Jo Jo 

Thus we find the solution as 

Aj>)-A« 



9(cos0') 
dcosfl' 

e(cos^) 



^ : #min — * 

6» : 7T - mi 



7T 



# } (35) 



where 



A<*>(0) 



A<>) 



-AW(0)+A,/2-A« 
A^(0)+A,-A« 

1 

F (arcsin?/ e , k e ) , 



g(l) 



"mm 



TT - 9 ■ 

in ~ * ) 



£-z\Jt-O z + 



(36) 



(37) 



and A« = Ai>S 1} ). 

Inverting ( l36l) . we derive cos 0(A) as 

cos6 l (A) = y / zTsn(w 6 i(A), ke), 



(38) 
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where 

[ 4K(ko)(\® + A 6 1] )/A e (0 < A (e) < A e /A - A«), 

u*(A) = 4K(k e )(-\M + A e /2 - A e l) )/A e (A fl /4 - A« < A^ < 3A./4 - A«), 

1 4K{ke)(\W ~A e + A^)/A e (3A„/4 - A« < A^ < A,). 

Using the results of (1351) . we derive t^ and 0^ as 



t (0) 



a 2 8z + 
1 



71 



E(k e )-E e (^ e ,k e ) 



2T e A (e) 

II fl (^ fl , ~Z_, kg) U(-Z-,kg) 



7T 



where ^ = arcsin[sn(M 9 , and 

E e (ij e ,k e ) — E(ip g , kg) — E(ip e (0), kg) 

for < A (6>) < A e /4 - A^\ 

= -Etye, h) + 2E(k e ) - S(Vto(0), k e ) 

for A /4 - A^ < A (e) < 3A e /4 - A^, 
= E(ij e , k e ) + A-E(kg) - £(<M0), M 

for 3A e /4 - A e 1] < A (e) < A e , 

U e (ip e , c, = Ufa, c, fe ) - n(^fl(0), c, 

for < A (9) < A fl /4- A^, 
= -U(ip e , c, fc e ) + 2n(c, k e ) - n(^(0), c, fc„) 

for A e /4 - A^ < A (e) < 3A e /4 - A { g\ 

= n(v>0, c, fe e ) + 4n( c , fc e ) - n(^(o), c, k 9 ) 

for 3A e /4 - A? } < A (9) < A fl . 



(39) 



(40) 



(41) 



4-2.2. dcos#(0)/dA < case In this subsection, we consider the case that the initial 
value of cos 0(A) satisfies dcos0(O)/dA < 0. We set 0(A = 0) = 9 K > in this subsection. 
Then A(0) in (154]) can be expressed as 
rcose dcos0' 

J cos oi 



g(2) 



e(cos^) 



+ 





cos £> 



+2 



o 



cos 


cos 9„ 



,(2) 



dcos0' 

e(cos^) 



+ 







cos 8. 



(2) 



+4 



COS £> n 



+ 



COS 



(2) 



10 JO 

Thus we find the solution as 

-A<»+Af 
\^(9) = { \V\e)+A 6 /2 + A e V 
-\V\9)+A e + A 9 2) 



dcos0' 
e(cos0') 
dcos0' 



e(cos0') 



3 (2) 



7T 



f : ( o 2) 



7T 



0„ 



0min 7 



9 ■ 9 min — > 9 



(2) 



(42) 



7T 



0„ 



mm "mm; 

■ ^0 {2) 
"mm ' "0 ' 



(43) 
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where \f ] = A^(^ 2) ). 

Then we obtain cos#(A) in the same form in ( l38i) inverting ( l4~3l) . but ug{\) in ( |38j) 
is modified as 

[ 4K(^)(-AW + AS 2) )/A 9 (0 < AW < A„/4 + Af), 

u e {\) = 4K(k ){\W - A e /2 - AS 2) )/A, (A fl /4 + A< 2) < A^ < 3A e /4 + A { e 2) ), 

{ 4K(ko){-\® +A e + Af Y )/Ae (3A e /4 + A?> < A< fl > < A e ). 

Using the results of (1421) . we derive t^ and as in the same form in (|39|) . but 

Eg(ipg, kg) and Ufl^, c, ke) are modified as 

^(^, k e ) = -E(4> e , kg) + £?(^(0), fc fl ) 

for < A w < A e /4 + A^ 2) , 
= + 2£(£; e ) + E(ip 9 (0), kg) 

for A,/4 + Af ] < A (e) < 3A„/4 + kf , 
= -£(^, k e ) + 4E{k e ) + EtyoiP), k e ) 

for 3A e /4 + Af ] < A (e) < A e , (44) 

n e (^e, c, = -n(^e, c, + n(^e(0), c, kg) 

for < A (e) < A e /4 + Af\ 

= n(^, c, + 2n(c, fc fl ) + n(^(o), c, 

for A e /4 + Af < A (e) < 3A e /4 + kf , 
= -n(^, c, + 4n(c, + 11(^(0), c, kg) 

for 3A e /4 + A^ 2) < A (e) < A e . (45) 

4-3. Consistency check of the analytical solution 

In this subsection, we compare the analytical results of bound geodesies, t(X), r(A), 
cos#(A) and <ft(\), with earlier literature as a consistency check. We can compare 



r(A) in (1261) with that in [T3J [22], in which the integral of motion, J dr/jR(r 



fdcos9/y 0(cos#), is solved in terms of Jacobi's elliptic function using # as the 
independent variable. Since [13] deals with null geodesies and [22] does not give explicit 
expressions including the turning points, we can not compare them exactly. However, 
we find that the formal expressions of r(A) in this paper and [T3], [22] are consistent. 
We can also compare the formal expression of cos6 l (A) in (J38l) with that in [13] . in 
which the integral of motion is solved in terms of Jacobi's elliptic function using r 
as the independent variable. Then we find that the formal expressions of cos#(A) in 
this paper and [13j are consistent. Although [13] derived t and <fi in terms of Carlson 
elliptic integrals [15] ; it seems difficult to compare the expressions of both t and in 
this paper with that of null geodesies in [IB] . Thus we compare the expressions of each 
integrals such as J* (r') 2 dr' / J R(r') in Appendix A with that in [15], in which formulas 



f° r ly T\i=i( a i + bit) Pl / 2 dt are derived in terms of Carlson elliptic integrals, where all 
quantities are real, x > y and ai + brf > for y < t < x. Using transformation of Carlson 



elliptic integrals of the third kind [23], it is straightforward to check P ' {r') 2 dr' / J R{r' 



in Appendix A agrees with the corresponding formula in [13]. We can also check the 
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other integrals in this paper agree with that in [J5]. These analytical checks show that 
the expressions in this paper are consistent. 

Moreover, we compare the analytical results in this paper with that of numerical 
integration method. It is a good check on the results that both the coefficients and 
turning points in each integrals are consistent. In numerical integration of geodesic 
equation, as explained in section [21 we transform r as r = pMj (1 + e cos ip) and cos 9 as 
cos 9 = cos 6^ cos x respectively. Using (jSj), we obtain the following set of differential 
equations [9] 



_ MyJ{\ -£ 2 ){(p-p 3 ) - e(p + p 3 cos ip)}{(p - p A ) + e(p -p 4 cos^)} 
dA = 1^2 • 

dx 



— = \]a?(l - £ 2 )(z + - z- cos 2 x), 

^ = T r (r) + T fl (cos0) + aC z , 
dA 

^ = $ r (r) + $,(cos 9) - aS } (46) 
dA 

where P3 = r%(l — e)/M andp 4 = r^(l + e)/M. We can numerically solve (1461) accurately 
without taking account of the turning points of both r and cos 9 because both ip and x 
are monotonic increasing functions of time. 

In figure [1], we compare the results of our analytical expressions with the results from 
the numerical integration method. We choose orbital parameters as a = 0.9M, p = AM, 
e = 0.7 and 9 inc = 40°. And set the initial values of ip and \ as ip(0) = and x(0) = 
respectively, which correspond to r(0) = = pM/(l + e) and 9(0) = 9^ = n/2. 
For numerical integration of ( 14*61) . we use the 4th order Runge-Kutta method with non- 
adaptive step-size control [T7] . This figure shows that the analytical solutions of geodesic 
equation in this paper exactly represent the solutions of bound geodesic orbits around 
a Kerr black hole. 



5. Summary 

We have derived analytical solutions for bound timelike geodesies in Kerr spacetime. 
This is the first time that analytical expressions of the fundamental frequencies are 
derived in terms of the elliptic integrals. The analytical expressions of the orbits, 
(t,r,cos9,<f)), have been also derived in terms of the elliptic integrals using Mino time 
as the independent variable for the first time. Since Mino time decouples the r and 9- 
motion, it leads to forms simpler than that in [13] for null case if we suitably transform 
variables, r and 9. We checked the consistency of the analytical expressions comparing 
them with the analytical expressions for the other cases, post-Newtonian approximation 
and numerical integration method. 

We can apply these solutions to the computation of gravitational waves from 
EMRIs. Gravitational waves from EMRIs are described by the Teukolsky formalism |24j. 
In the frequency domain calculation of the Teukolsky formalism [25l [26] , we can use 
the analytical solutions directory [271 EH] and compute the orbits more accurately than 
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500 1000 1500 2000 2500 



t/M 

Figure 1. Comparison of the function x(t) = r(t) sin 0(f) cos 4>{t) computed using the 
analytical expressions of this paper with the result of a numerical integration. In this 
figure, we set orbital elements as a = 0.9M, p = AM, e = 0.7 and 9 mc = 40°. And we 
set initial value of ip and x as V'(O) = and x(0) = respectively, which correspond 
to r(0) = r { 1] = pM/(l + e) and 9(0) = 9^ = tt/2. Upper figure shows plots of both 
analytical solution, XA(t), and numerical integration method, x^(t). Lower figure 
shows the residual errors between the results of both analytical solution and numerical 
integration method, XA(t) — x^(t). 

numerical integration of geodesic equation. In principle, we can compute the orbits with 
machine accuracy. Using the analytical expressions of radial and polar motion in this 
paper, we showed that the analytical expressions enable one to compute gravitational 
waves from EMRIs very accurately [28]. Although it may takes longer time to compute 
the orbits using the analytical solutions than using numerical integration method [12J, 
but see [13], it is not serious in computing gravitational waves. This is because we 
compute orbits only for one orbital period of radial and polar motion, A r and Ag, and 
computation time of the orbits, ~seconds, is much smaller than that of gravitational 
waves from EMRIs, ~ hours to days [271 [281 29J. Thus we believe that the analytical 
solutions are very useful for the computation of gravitational waves from EMRIs. In 
the time domain calculation of the Teukolsky formalism (see brief review in section 
3.8 in [10]), we may need the inversion of t(X) in order to compute the orbits, t(A), 
cos#(A), t(X) and 0(A), in the coordinate time. Although we do not know the analytical 
expression of A(t), we may easily obtain A(t) by numerical iteration if we set the initial 
solution as A = t/T. Thus it may also be useful in the time domain calculation if 
the numerical iteration converges faster than the numerical integration of the geodesic 
equation. 

We may also apply these solutions to investigate the properties of geodesies of Kerr 
black holes. Although it seems difficult to classify orbits in the strong field because 
of its complexities, Levin et al recently suggested a taxonomy of orbits introducing 
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a rational number which is constructed from orbital frequencies [SUl EE]- Both the 
analytical expressions of the fundamental frequencies in this paper and the taxonomy 
of orbits may help us to discuss the conditions characterizing zoom-whirl orbits [6] and 
other extreme phenomena in Kerr backgrounds. The other applications may be null or 
unbound geodesies. We can apply our method to them with a few modifications. For null 
geodesies, we have to eliminate the mass term of the small body in the geodesies. For 
unbound geodesies, we can express the orbits in terms of the elliptic integrals using Mino 
time although there are no fundamental frequencies for the orbits. However, it should be 
noted that we may have to improve computation time when we consider null or unbound 
geodesies because we have to trace the orbits for longer time than bound orbits cases. 
If we can not improve computation time, we may have to use both analytical solutions 
and numerical integration [12]. Finally, we note that we can not use the analytical 
solutions in this paper when the inclination angle from the equatorial plane of black 
hole is 6> inc = ir/2. This is because &q(cos9) in ([I]) diverges when 9 = 0, and its elliptic 
integral also diverges. We do not know how to address this issue without any further 
approximation though one can solve it if one uses a post-Newtonian expansion |21j. All 
of them will be discussed in a future work. 
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Appendix A. formulas of integrals of radial motion 

In this appendix, we derive some formulas which are needed to obtain T t (r) and T e w in 
section Eland t^ r ' and <fr- r > in section HI In order to compute them, we have to investigate 
/; dr'/{(r' - r±)^R(P)}, f^r' fdr' / ^R^) and g dr'/{(r' - M) 2 /ff(7)}, see (119]). 

Since r'dr' / J ' R{r') is derived in sectionEland J r r 2 dr'/ { (r'—M) J R(r') } can be derived 

from / r r 2 dr'/{(r' — r±)y i?(r')} when we set r± = M, i.e. a = M, we do not show again 
these expressions in this appendix. 

As we derived r'dr' / J R(r') in section El it is useful to transform r into y r . Then 
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we have the following relations. 



1 



r 2 - r 3 



1 



r 2 - r± 1 - h±y 2 r 



r — r± r 3 — r± 

2 2 , ( r 2 - r 3 )(r 2 + 3r 3 ) 1 (r 2 -r 3 ) 5 
t = r 3 H ~ ~ ; — n + 



t-3- M y 

r - M ) 



1 + ^ 



2 1 - Kyi 

1 r 2 - r 3 ( r 2 -r 3 
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(y r - hr 1,2 Y (y r + hr 1/2 f\ 



2r 2 -M\r 2 -M J 1 - h M y 2 
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(A.l) 
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r 2 - r 3 



4:h M \r 2 - M 



1 



where fr M = h±(a = M) = (n - r 2 )(r 3 - M)/[(ri - r 3 )(r 2 - M)]. 

Then it is straightforward to compute J r r 2 dr'/{(r' — r±)^ii(r')} as 

dr' 2 



r 2 ( r ' _ r± )^R(r') (r 3 - r±) ^/(l - £ 2 ) (r x - r 3 ) (r 2 - r 4 ) 

r 2 - ^3 



[F(arcsiny r , k r 



-n(arcsiny r , —h±, k r 



(A.2) 



r 2 - r± 

However, we need a reformulation of the last terms of both / r r 2 (r') 2 dr' / y R(r') 

and /; 2 dr'/{(r' - M) 2 /r(7)}. If we set J n [c] = f%dy'/{{y' - cfy^/)}, where 
= (1 — y 2 ){l — k 2 y 2 ), we can represent these terms as J 2 [c] + J 2 [— c]. Using reduction 
formula of the elliptic integrals, see section 17.1.5 in [14], we can express J 2 [c] + J 2 [— c] 
in terms of the elliptic integrals as 

Mc] + J 2 [-c] = — 7-7 { [(2c 2 - l)/c 2 - l] n(^, -c" 2 , Av) + (1 - k*c 2 )F{4>, k r ) - Ety, k r ) 



<f(c) 



y 2 — c 2 



(A.3) 



where ip = arcsiny. 



Then we can express JI (r') 2 dr' / J R(r') and f* ' dr'/{(r' — M) 2 J R{r')} respectively 



as 
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Appendix B. \a\ = M case 

In this appendix, we show the analytical expressions of T, *W and in the case 



\a\ = M. Using the results of Appendix A , we derive them as 

2a 2 £ z + T t 
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